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In the prototype sandpile model of self-organized criticality 
time series obtained by decomposing avalanches into waves of 
toppling show intermittent fluctuations. The q-th moments 
of wave size differences possess local multiscaling and global 
simple scaling regimes analogous to those holding for velocity 
structure functions in fluid turbulence. The correspondence 
involves identity of a basic scaling relation and of the form 
of relevant probability distributions. The sandpile provides a 
qualitative analog of many features of turbulent phenomena. 
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Intermittency is a characterizing feature of turbu- 
lence • Intervals of low activity and of irregular du- 
ration are separated by sudden bursts in a wide range 
of strengths. Random quantities like activity variations 
across fixed space or time distances, are distributed ac- 
cording to non-Gaussian laws, whose moments display 
also a peculiar multiscaling. At first sight, some qualita- 
tive features of turbulent signals seem to be recognizable 
in the distant context of self-organized critical dynam- 
ics ||. In this field suitable models describe the bursts 
of activity by which a system, subjected to an external 
inflow of energy or matter, relaxes under the control of 
local, nonlinear threshold mechanisms. The bursts are 
called avalanches, and their time series show indeed some 
similarities to those of turbulent fluctuating quantities. 
Of course, these similarities are only apparent. Indeed, 
as a rule, different avalanches are just uncorrelated, and 
the probability distribution functions (pdf 's) of avalanche 
quantities do not show multiscaling. 

Quite recently it was realized that the two-dimensional 
(2D) Bak-Tang-Wiesenfeld (BTW) sandpile, the proto- 
type model of self organized criticality , obeys a pecu- 
liar form of multiscaling for the pdf's of several avalanche 
measures §. This result came after many years of work 
assuming allmost unanimously the simple finite size scal- 
ing form. The unexpected multiscaling makes the 2D 
BTW one special among models in the same class, and 
raises the issue of understanding its origin and the pos- 
sible analogies with other complex scaling phenomena. 

In this Letter we show that one can indeed establish a 
surprisingly close correspondence between BTW sandpile 
dynamics and fully developed turbulence. This corre- 
spondence includes an identity of the basic relations con- 
necting the multiscaling exponents of the relevant pdf's. 

The BTW model is denned on a square lattice (L x 
L) box Q. On site i z% = 0, 1, 2, . . . is the number of 



"grains". If z, < 4, Vi, the configuration is stable. Grain 
addition to a stable configuration is made by selecting 
at random a site k where Zk — > Zk + 1. If then Zk > 4, 
toppling occurs, i.e. z% — > Zk — 4, while each nearest 
neighbor, I, of site k gets one grain (zi — > zi + 1). If A; 
is at the border grains are dissipated. The toppling of 
site k may cause instabilities in the neighbors, leading to 
further topplings at the next microscopic time step, and 
so on. Thus, an avalanche made by a total number s > 
of topplings occurs before a new stable configuration is 
reached and a new grain is added. After many additions 
the system reaches a stationary critical state in which 
avalanche properties are sampled. 

A key notion in the approach to the BTW is that of 
wave decomposition of avalanches The first wave is 
obtained as the set of all topplings which can take place 
as long as the site of addition is prevented from a possible 
second toppling. The second wave is constituted by the 
topplings occurring after the second toppling of the addi- 
tion site takes place and before a third one is allowed, and 
so on. The total number of topplings in an avalanche is 
the sum of those of all its waves. The global sample of all 
waves resulting from the decomposition of a sequence of 
several avalanches has exactly known scaling properties, 
well confirmed numerically |6|,|^] . The wave size has a pdf 
satisfying finite size scaling P w (s,L) ~ s _T ™ f w (s/L D "'), 
where, in d = 2, t w = 1, D w = 2, and f w is a suitable 
scaling function. Most recently it has been shown that 
the time series of successive waves sizes {s m },m = 1,2,.., 
has long range autocorrelation in the 2D BTW sand- 
pile 0. This is not the case for other similar models, like 
the Manna one ||, in which the size of a wave is just un- 
correlated with that of the preceeding one [|| . The long 
time autocorrelation of waves is responsible for the fact 
that the scaling of the BTW avalanche size pdf, P av , is 
different from that of P w ||. It remains a major chal- 
lenge to identify which features of the BTW wave time 
series are at the root of the multifractal scaling pattern 
detected for avalanches ||. Here we find a key to this 
issue by establishing a correspondence between g-th or- 
der wave size difference moments of the sandpile, and 
velocity structure functions of a turbulent fluid flow. 

In fully developed turbulence one measures a ve- 
locity field u(x), with structure functions G q (r) = 
(\u(x + r) — uix)^) 1 ^ G q displays a power law r- 

dependence G q ~ r^ q in the inertial range 1 < r/rj < 
£ ~ (l/ v ) ~ i? 3 / 4 , where £ is a correlation length, R is 
the Reynolds number, and I and r] are the integral and 
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the dissipation scales, respectively. The local dissipation 
rate e ~ (du/dx) 2 has anomalous fluctuations. Indeed, 
the values of du/dx have a broad pdf, with moments 
obeying multiscaling in £. By matching the scaling of 
\du/dx\ with that of \u(x+rj) — u(x) \ in the inertial range, 
one also expects G q (r ~ rf) ~ £ a<! , where a g is another q- 
dependent exponent |yj . On the other hand, as one con- 
siders u-increments for r/r) 3> £, the corresponding pdf 's, 
due to the lack of correlation, become r-independent and 
give information only on the global sample of velocity 
variations over all r, i.e. G q (r) ~ (\u — (u) \ q ) 1 ^ q ', for 
r/rj^^. The pdf's of u(x + r) —u(x) become also Gaus- 
sian while multiscaling crosses over to a global simple 
scaling, so that G q (r) ~ £ Q , for r/rj 3> £, independent of 

The scaling of G g is often detected by measuring the 
velocity field u at different times at a fixed position in 
space Q. The data of this time series are then con- 
verted by interpreting time intervals as space intervals, 
once the average convection velocities are known. This 
suggests to regard the time series {sk} of BTW wave 
sizes in a similar perspective: also in this case we deal 
with a long range correlated signal, and it is legitimate 
to ask whether turbulence-like mechanisms in this signal 
are responsible for some form of multiscaling. Of course, 
a conversion of time into space dependence would not 
make sense in this case. Thus, in the analogy t will def- 
initely replace r. A structure function for the sandpilc 
can be defined as 



F q (t) = (\s k+t -s k \ q ) 1/q , 



(1) 



In turbulence the difference in velocities is a measure of 
how a perturbation of the average flow grows by prop- 
agating from one point to another. In sandpile dynam- 
ics the difference in size between two successive waves 
can also be considered as a measure of how a pertur- 
bation evolves in time. The instantaneous entity of the 
perturbation is represented by the number of topplings 
occurring in the corresponding wave. In analogy with 
turbulence, a scaling form one can expect for F q is 



F q (t) ~ C"t^f g (t/o, 



(2) 



where f q are scaling functions, and a q and Cq are expo- 
nents possibly varying with q, and specified here by the 
same symbols as the corresponding turbulence ones. £ 
in the BTW case is a correlation time for waves, which 
should scale ~ L z . It plays the role of the correlation 
length in turbulence, and is denoted by the same symbol 
here. The long time correlation is a peculiar feature of 
2D BTW wave series, essential for establishing the cor- 
respondence with turbulence. Avalanche time series in 
general show finite correlation time for models of self or- 
ganized criticality, including the 2D BTW sandpile |To| ]. 

If assumed for G q , the form (2) embodies the various 
scaling regimes described above for turbulence fLN. la- 



the t dependence when i ^ £ is correctly absorbed. In 
this limit one gets then F q ~ ^ a i+Cq , We know further 
that, for t ^> £, F q should reproduce the moment scaling 
of the global wave distribution P w . Indeed, wave sizes 
are totally uncorrelated for t 3> £, so that F q [t 3> £) ~ 

{\s-{s) P J q fl\ with {-) Pw ~ J-P w (s,L)ds. Since 

P w has a simple, constant gap (D w = 2) scaling form, 
F q (t » £) ~ L 2 . Thus, the exponents in Eq. (g) must 
satisfy the identity 



a q + ( q = 2 /z 



(3) 



for all q. A scaling relation analogous to 3 is expected 
to hold for the corresponding turbulence exponents [ pT| . 
For 1 < t < £ one expects f q ~ const., which leads to 
a power law i-dependence, F q (t) ~ t^ q , in the analog of 
the inertial range. Finally, for t ~ 1, which corresponds 
to r ~ r), one recovers the moments of the equivalent 
of \du/dx\ in turbulence: F q (l) ~ . The interval 
between successive waves, taken here as the unit of time, 
corresponds to rj in the turbulent case. 

By analyzing samples of up to 10 7 waves for L = 
128, 256, 512 and 1024, we could verify that the scaling 
ansatz 2 and relation 3 are very well satisfied. Fig. [j] re- 
ports collapse plots of the quantity y^ q f q {y) for various 
g's. The collapse quality increases with q and z ~ 0.7 is 
an optimal choice for all cases. Since f q is approximately 
constant for y < 1, from the low y part of the plots one 
also obtains estimates of £ g . The i-dependence of F q in 
the inertial range shows some degree of extended self- 
similarity, a well known property of structure functions 
in turbulence @ : log-log plots of F q (t) q vs Fx (t) (Fig. g) 
are rather straight in the range corresponding to the in- 
ertial regime. In spite of the relatively narrow interval of 
variation of the moments in this range, this allows rather 
precise determinations of Cg/Ci- These, combined with 
an estimate of £i from the plots in Fig. |l| lead to the 
values reported in Table Q. The exponents a q reported 
there are determined by log- log plots of F q (l) versus L, 
which arc definitely linear and show clear g-dependent 
slopes. The values of z( q and za q reported in Table || 
satisfy Eq. (^) very well. 

A key fluctuating quantity in the description of tur- 
bulent flow is the dissipation rate averaged over spatial 
volumes of linear size r, e r [0 . In order to get further in- 
sight into the correspondence with turbulence, we studied 
an analogous wave quantity, given by the average over an 
interval t of what should correspond to the local dissipa- 
tion rate: 



t 

E 

fc=i 



(sk+i - s k ) 2 /t 



(4) 



deed, we expect f q (y — > oo) 



y 



-c 9 



in such a way that 



In analogy with turbulence we expect (ef) / (e t ) q ~ t~^ q , 
in the inertial range. Fig 3 shows plots of this quantity 
which clearly manifest multiscaling and are remarkably 
similar to those one obtains for e r in turbulence [2|. The 
plots of course do not include the plateaus one finds in the 
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turbulent dissipation range (r < 77), which is not covered 
by our correspondence. Convergence to zero for all q 
can also be perceived for large t. These results further 
support the identification of Sk+i — Sk as the analog of a 
velocity gradient in turbulence. 

Thus, F q in the BTW sandpile has a scaling isomor- 
phic to that of G q in turbulence. So far, theoretical 
approaches always tried to link avalanche to wave scal- 
ing |l^,[u|. Since this last scaling is simple and ex- 
actly known at global level, attempts concentrated on the 
connection between avalanche scaling and this regime, 
just overlooking the very possibility of a different local 
regime and its consequences. We find that, in the range 
1 < t < L z , the wave series displays complex multiscal- 
ing structure, of the same type as that found for velocity 
fields in turbulence. This is not the case for the Manna 
sandpile || , which does not show any form of multiscal- 
ing. In that model, due to the absence of correlation, the 
global regime of the structure function is reached as soon 
as t = 1 and there is no window for local and inertial 
multiscaling regimes. 

The intermittent local wave regime is certainly at 
the basis of the multiscaling already detected for global 
avalanche distributions of the 2D BTW model 0. An 
open problem remains that of connecting the multifractal 
exponents of the wave time series with those of avalanche 
distributions. Wave and avalanche time series are sepa- 
rated by an infinite coarse graining in time, for L — > 00. 
At numerical level, one could in principle perform this 
coarse graining by partitioning waves into successive sets 
of n elements Q. The corresponding block variables 
Sk = s n (fc-i)+i + s n (k-i)+2 + ■ • • + s n k, k = 1, 2, . . . are 
expected to have scaling properties approaching those 
of avalanche sizes in the limit n — > 00. So, for coarse 
grained structure functions F q , n = (\Sk+t — Sfcl 9 ) 1 ^ 9 the 
scaling of Eq. (|2|) should cross over, for increasing n, to 
the t- independent form F q ^ n ~ L aq / q , where a q is the 
multifractal moment exponent of the global avalanche 
size distribution P av (s, L), i.e. {s q ) Pav ~ L aq ||. Unfor- 
tunately, sampling limitations make such coarse graining 
extremely difficult to realize for large enough n. 

The analogy with turbulence goes beyond the scaling 
properties alone. The pdf 's of velocity differences in tur- 
bulence show an increasing asymmetry between positive 
and negative increments as r is reduced within the iner- 
tial range 0. This asymmetry, implying nonzero skew- 
ness, is an important feature of energy transfer across dif- 
ferent scales. Fig. || reports a plot of the pdf of Sk+i — Sk 
for L = 1024, which shows the expected slight asym- 
metry. One can also try to apply to such plots fitting 
schemes used in turbulence. We adopted a scheme pro- 
posed most recently, based on non-extensive statistics as- 
sumptions, and extremely successfull in reproducing fluid 
turbulence experimental data p5| . Omitting details, to 
be published elsewhere, we just include in Fig. |] the ob- 
tained best fit. This is very satisfactory within a sub- 
stantial deviation range and the two parameters involved 



obey within less than a percent a relation established in 
Ref. ||. 

Summarizing, BTW wave time series produce multi- 
scaling schemes and pdf's isomorphic to those of fully 
developed turbulence. This identifies the origin of the 
peculiar scaling of the 2D BTW sandpile, compared to 
similar models, and gives it a novel theoretical valence. 
Further elucidating why BTW waves behave like veloci- 
ties of a turbulent flow, or why (sk+i — Sk) 2 can be viewed 
as a fluid energy dissipation rate, should provide deeper 
insight into this problematic model. Multiscalings con- 
sistent with Eqs. (2-3) have been previously identified 
also for a model of epitaxial growth in ID |ll|] and, most 
recently, for the directed polymer in tilted columnar dis- 
order in 2D |u|. In the second case careful determina- 
tions of a q and £ q revealed that the multiscaling is in fact 
just a bi-scaling, of the kind expected for Burgers turbu- 
lence |T^], rather than for turbulence plpl. Our results 
seem to definitely exclude bi-scaling in favor of genuine 
turbulent multiscaling. Thus, the analogies with turbu- 
lence discovered here include this crucial aspect. A pos- 
sible connection between turbulence and self organized 
criticality was first conjectured when this paradigm was 
introduced 0] . Our results show that indeed there exists 
a precise analogy in the case of the BTW sandpile, and 
that this analogy can be established through the concept 
of waves, which are in turn peculiar to a restricted class 
of models. However, this does not exclude that similar 
links with turbulence could be established on a different 
basis fl8f| or also for other models, in which for example 
waves are not defined p9[ . 

Recent studies of the statistics of solar flares have 
correctly emphasized substantial differences existing be- 
tween self-organized criticality and turbulence mecha- 
nisms, for the explanation of power law pdf's pC| . Such 
differences, which favor turbulence models, can ulti- 
mately be ascribed to the lack of long time correlations 
for the avalanches in the self organized critical models 
used so far to describe these phenomena. Quite re- 
markably, the distinction between the two mechanisms 
is however lost in the 2D BTW model, which has long 
time correlated bursts and genuine turbulence, if stud- 
ied at the wave time scale. Self-organized critical models 
in which activity bursts (no matter whether waves or 
avalanches [Tc|]) have long time correlations, would cer- 
tainly be better candidates to describe flares. Of course, 
these correlations alone would not be necessarily enough 
to produce also turbulent scaling patterns. 

Partial support from the European Network No. 
ERBFMRXCT980I83 is acknowledged. 



[1] U. Frisch, "Turbulence" , (Cambridge University Press, 
Cambridge, 1995). 



3 



[2] D. Lohse and S. Grossmann, Physica A 194, 519 (1993). 
[3] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 

59, 381 (1987); Phys. Rev. A 38, 364 (1988). 
[4] D. Dhar, Physica A 263, 4 (1999). 

[5] C. Tebaldi, M. De Menech, and A. L. Stella, Phys. Rev. 

Lett. 83, 3952 (1999); M. De Menech, A. L. Stella, and 

C. Tebaldi, Phys. Rev. E 58, R2677 (1998). 
[6] E. V. Ivashkevich, D. V. Ktitarev, and V. B. Priezzhev, 

Physica A 209, 347 (1994). 
[7] D. V. Ktitarev, S. Liibeck, P. Grassberger, and V. B. 

Priezzhev, Phys. Rev. E 61, 81 (2000). 
[8] M. De Menech and A. L. Stella, Phys. Rev. E 62, R4528 

(2000). 

[9] S. S. Manna, J. Phys. A 24, L363 (1991). 
[10] For an exception se e, e.g., J. Davidsen and M. Paczuski, 

|cond-mat/0105532 . 
[11] J. Krug, Phys. Rev. Lett. 72, 2907 (1994). 
[12] V. B. Priezzhev, E. V. Ivashkevich, and D. V. Ktitarev, 

Phys. Rev. Lett. 76, 2093 (1996). 
[13] M. Paczuzki and S. Boettcher, Phys. Rev. E 56, R3745 
(1997). 

[14] R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, and S. 
Succi, Phys. Rev. E 48, R29 (1993). 

[15] C. Beck, G. S. Lewis, and H. L. Swinney, Phys. Rev. E 
63, 035303(R) (2001). For non-extensive statistical me- 
chanics, see C. Tsallis, Braz. J. Phys. 29, 1 (1999). 

[16] R. Mohayaee, A. L. Stella, and C. Vanderzande, Phys. 
Rev. Letters, 87, 085701 (2001). 

[17] See, for example, J. P. Bouchaud, M. Mezard, and G. 
Parisi, Phys. Rev. E 52, 3656 (1995). Bi-scaling would 
imply, e.g., a q(, q piece-wise linear in q, with just two 
distinct slopes. 

[18] A. Erzan and S. Sinha, Phys. Rev. Lett. 66, 2750 (1991). 
[19] The forest fire model [P. Bak, K. Chen, and C. Tang, 

Phys. Lett. 147, 297 (1990)] has been recently indicated 

as a good candidate: K. Chen and P. Bak, Phys. Rev. E, 

62, 1613 (2000). 
[20] G. Boffetta, V. Carbone, P. Giuliani, P. Veltri, and A. 

Vulpiani, Phys. Rev. Lett. 83, 4662 (1999). 




-0.1 1 ■ 1 ■ 1 1 

4.15 4.35 4.55 

log 10 F,(t) 

FIG. 2. Extended mutual power law scaling between dif- 
ferent moments of the pdf of \s^+t — Sk\- 
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FIG. 3. Log-log plots (e 9 ) / (et) 9 and estimates of fi q . 




FIG 1 Collapse plots for f FIG. 4. Distribution of the nomalized wave size difference 

(solid line), and relative fit (empty circles). 
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ZOt q 


0.5 


0.46 


1.43 


0.8 


0.37 


1.58 


1.0 


0.33 


1.64 


2.0 


0.21 


1.79 


3.0 


0.16 


1.85 


4.0 


0.13 


1.88 


6.0 


0.091 


1.91 



TABLE I. For each q the uncertainties in both quantities 
are of the order of the discrepancy with respect to Eq. []. 



